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CONSTRAINING PROPOSED COMBINATIONS OF ICE HISTORY AND EARTH RHEOLOGY 
USING VLBI DETERMINED BASELINE LENGTH RATES IN NORTH AMERICA 

J. X. Mitrovica 1 , J. L. Davis 2 , and I. I. Shapiro 2 


Abstract . We predict the present-day rates of change 
of the lengths of 19 North American baselines due to 
the glacial isostatic adjustment process. Contrary to 
previously published research, we find that the three- 
dimensional motion of each of the sites defining a base- 
line, rather than only the radial motions of these sites, 
needs to be considered to obtain an accurate estimate 
of the rate of change of the baseline length. Predic- 
tions are generated using a suite of Earth models and 
late Pleistocene ice histories; these include specific com- 
binations of the two which have been proposed in the 
literature as satisfying a variety of rebound related geo- 
physical observations from the North American region. 
A number of these published models are shown to pre- 
dict rates which differ significantly from the VLBI ob- 
servations. 

Introduction 

Glacial isostatic disequilibrium is manifested in a num- 
ber of geophysical observables, none more apparent than 
the associated three-dimensional crustal deformations. 
Radial deformations have received the most attention 

I e.g. Cathles, 1975]. Horizontal motions, in contrast, 
lave been considered in only a very few studies [e.g., 
James and Morgan , 1990; Mitrovica et aL, 1993a], since 
there has, until very recently, been an absence ot obser- 
vational constraints on such motions. The continuing 
improvement of space-geodetic techniques has enabled 
three-dimensional crustal deformation rates to be esti- 
mated with an accuracy useful for application to the 
glacial isostatic adjustment problem. Indeed, Mitro- 
vica et al. [1993a] compared theoretical predictions of 
three-dimensional motions associated with the adjust- 
ment process with available very-long- baseline interfer- 
ometry (VLBI) constraints from North America and 
Europe. In particular, Mitrovica et al. [1993a] used 
the results of Ryan et al. [1993] who analyzed an en- 
semble of Mark HI geodetic VLBI data obtained from 
1979-1991 to determine the evolutions of a large num- 
ber of baselines. The “global” solution of Ryan et al 
[1993], which they termed GLB868, included 864,359 
group delay measurements obtained in 1648 observing 
sessions at a globally distributed network of sites. 

The typical application of observables related to gla- 
cial isostatic adjustment is to constrain the radial pro- 
file of mantle viscosity. A complication in most such 
applications is the sensitivity of the geophysical data 
to the space and time geometry of the late Pleistocene 
ice sheets. Mitrovica et al. [1993a] concluded that the 
VLBI constraints are most usefully applied to assess the 
acceptability of specific ice model /Earth model pair- 
ings. In this paper our main goal is to use the con- 
straints on baseline velocities provided by the GLB868 
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solution to test the acceptability of a set of published 
models for the combination of ice history and viscosity 
profile associated with the North American region. 

Method and Results 

In the present analysis we focus on the rates of change 
of bn line lengths (hereafter “length rates”) for sites 
in Ni rth America. We limit our attention to base- 
lines which extend, at least in part, over the near field 
and periphery of the ancient Laurentide ice sheet. We 
will not use baselines which include sites on the west 
coasts of Canada and the contiguous United States, in 
order to avoid areas of well-known tectonic deforma- 
tions, nor baselines which include the Ft. Davis VLBI 
site, which are characterized by anomalous, and as yet 
unexplained, variations in length. Applying these crite- 
ria to the sites included in the GLB868 solution yields a 
set of 19 baselines; these are listed in Table 1, (together 
with the length rates obtained in the GLB868 solution) 
and illustrated in Figure 1. (Ryan et al. [1993] scaled 
the uncertainties so that the reduced x 3 for the series of 
baseline estimates is unity. For Richmond- Westford we 
have used two significant digits for the a because round- 
ing errors would be large with respect to the value.) 

Our predictions of baseline length rates will be based 
on the spectral formalism for computing three dimen- 
sional motions detailed in Mitrovica et al. [1993b], which 
is appropriate for the surface loading of a spherically 
symmetric, self-gravitating, linear visco- elastic Earth 
model. We adopt a spherical harmonic truncation level 
at degree and order 256 for the surface mass load. The 
load will be comprised of a model for the late Pleis- 
tocene ice sheets, and a gravitationally self-consistent 
ocean loading component. The latter is computed us- 
ing the pseudo-spectral method of Mitrovica and Peltier 
[19911 . The Earth models we consider all have an elas- 
tic structure given by the Earth model PREM and a 
Maxwell visco-elastic rheology. The set of radial viscos- 
ity profiles used in the calculations are specified below. 
In general these profiles are defined by the thickness of 
an elastic lithosphere (denoted LT), and by a constant 
value for each of the viscosities in the upper and lower 
mantle, respectively, v UM and v LM . (The boundary be- 
tween the upper and lower mantle region is taken to be 
at 670 km depth.) Our results will often be summarized 
in terms of the x 3 residual per degree-of-freedom (here- 
after “reduced x 3 ”) of the VLBI baseline length rates 
relative to the predictions obtained using the specified 
combination of ice history and radial viscosity profile. 
The reduced x 3 computed assuming that glacial iso- 
static adjustment does not influence the baseline length 
rates listed in Table 1 (hereafter “the reduced x 3 of the 
null hypothesis”) is 3.3. 

In the following discussion we refer to a “standard” 
Earth model which has LT=120 km, u UM = 10 21 Pas and 
v LM = 2 x 10 21 Pa s. In Figure 2 (solid line) we sum- 
marize predictions of the length rate of the Richmond- 
Westford baseline for a suite of Earth models identical 
to the standard model with the exception that u LM is 
varied from 5 x 10 20 Pa s to 5 x 10 21 Pa s. The ice loading 
used in the calculation was based on the recently pub- 
lished ICE-3G deglaciation chronology, which extends 
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TABLE 1. GLB868 Baseline Length Rates 


Baseline Rate (mm/yr) 


Algonquin Park (A)-Gilmore Ck (G) 
Algonquin Park (A)-Westford (W) 
Gilmore Ck (G)-Haystack (H) 
Gilmore Ck (G)-NRAO 85 (Nl) 
Gilmore Ck (G)-Platteville (P) 
Gilmore Ck (G)-Richmond (R) 
Gilmore Ck (G)-Westford (W) 
Gilmore Ck (G)-Whitehorse (Wh) 
Beltsville (B)-NRAO 85 (Nl) 
Beltsville (B)-Richmond (R) 
Beltsville (B)-Westford (W) 
Haystack (H)-NRAO 140 (N2) 
Maryland Pt (M)-Richmond (R) 
Maryland Pt (M)-Westford (W) 
NRAO 85 (Nl)-Richmond (R) 
NRAO 85 (Nl)- West ford (W) 

NRAO 140 (N2)-Westford (W) 
Platteville (P)-Westford (W) 
Richmond (R)-Westford (W) 


3.3 ±0.8 
-0.7 ±0.5 
-2.3 ±2.2 
-0.5 ± 1.4 

5.6 ±2.3 
1.1 ± 1.2 
-0.4 ±0.5 
-3.8 ±3.8 

3.4 ± 1.3 
6.9 ±2.5 
0.0 ±3.0 

1.4 ±2.0 
0.6 ±4.4 
1.4 ±1.5 

2.2 ±0.9 
4.0 ±2.0 
0.5 ±0.3 

3.3 ±2.2 
-0.2 ±0.15 


from 18 to 5 kyr before present ( Tushingham and Peltier 
[1991]), and adapted in order to incorporate a growth 
phase for the late Pleistocene ice sheets (see Mitrovica 
and Peltier [1993]). 

The analysis of Tushingham [1991] represents the only 
attempt to date to calculate baseline length rates due to 
glacial isostatic adjustment. However, he assumed that 
these rates could be accurately computed by consider- 
ing only radial deformations at the terminating sites. 
To assess the validity of this assumption, we also show 
in Figure 2 (dashed fine) the component of the baseline 
length rate which is due to radial motions at the Rich- 



Fig. 1. The VLBI baselines included in this study. The 
letters on the figure denote specific VLBI sites, as indi- 
cated in Table 1. The shaded region on the plot shows 
the geographic extent of the North American ice cover 
during the last glacial maximum according to the ICE- 
30 model of Tushingham and Peltier [1991]. 


Richmond - Westford 



Fig. 2. The predicted length rate of the Richmond- 
Westford baseline as a function of The solid line 

was computed from the three-dimensional motion of the 
Richmond and Westford sites, while the dashed line was 
generated assuming that only radial motions at these 
sites contribute to the length rate. The shaded region 
represents the VLBI determined constraint (±la) on the 
baseline length rate (Table 1). 


mond and Westford sites. It is clear that the Tushing- 
ham [1991] assumption is, in general, invalid. Indeed, 
over the range of Earth models considered in Figure 2 
the horizontal motions of the end sites generally con- 
tribute more to the baseline length rate than do the as- 
sociated radial motions. Clearly, inferences of viscosity 
based upon predictions generated using the Tushingham 
[1991] approximation would be suspect. 

In Figure 3 we plot the reduced y 2 for the comparison 
of the VLBI baseline length rates with the predictions 
generated from many Earth models and two ice mod- 
els. Figure 3a provides results generated using the ice 
loading adapted from the ICE-3G deglaciation chronol- 
ogy; each curve connects results from a set of Earth 
models constructed by varying a specific characteristic 
of the standard model (see caption). The two curves in 
Figure 3b are distinguished by the ice loading history 
adopted in the calculation. In particular, we generated 
the points on the dashed-dotted line using the ICE-1 
model of Peltier and Andrews [1976] (adapted in the 
manner described above). 

The results in Figure 3 show that the Earth models 
and icc histories represented there do not yield predic- 
tions with reduced y 2 values significantly less than the 
value for the null hypothesis (3.3). The reasons for this 
outcome may be some combination of errors in the de- 
tails of the ice loading histories over North America, 
the limited class of Earth models considered (see be- 
low), and contributions from neglected geophysical pro- 
cesses (e.g., tectonic deformations). Nevertheless, two 
important conclusions may be made. First, some ice 
history/Earth model parings predict length rates which 
appear to be incompatible with the VLBI determined 
constraints, and may as a consequence be ruled out. For 
example, variants of the standard Earth model with LT 
< 80 or v lM > 3.3 x 10 21 Pas, together with the ICE-3G 
dcglaciation chronology, yield reduced y 2 residuals in 
excess of 6.6 (or larger than twice the value for the null 
hypothesis). In contrast, many Earth models do not 
produce extremely high reduced y 2 values, and there- 
fore cannot be rejected. 
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Fig. 3. The reduced x 4 residual for the VLBI deter- 
mined baseline length rates in Table 1 relative to the 
corresponding predictions obtained from a suite of Earth 
models and ice histories, (a) All curves are computed 
using the ice loading history adapted from the ICE-3G 
deglaciation chronology. The solid, dotted, and dashed 
lines are generated by varying, respectively, v LM ,v VM , 
and LT of the standard Earth model. Each curve has 
an associated abscissa scale, (b) Both curves are gener- 
ated by varying u LM of the standard Earth model. The 
solid and dashed-dotted lines were computed using ice 
loading histories adapted from, respectively, the ICE-3G 
and ICE-1 deglaciation models. 


The predicted length rates for a subset of six of the 
19 baselines, for each family of results shown in Fig- 
ure 3a, are given in Figure 4. The dramatic increase 
in the reduced x 2 i for Earth models with progressively 
thinner lithospheres, is clearly the consequence of a mis- 
fit between predictions and observations for almost all 
the baselines considered in Figure 4, suggesting a rather 
systematic preference for Earth models with relatively 
thick lithospheres (i.e., LT > 100 km). The pronounced 
increase in the reduced x% for v LM greater or less than 
a value near 10 31 Pa s, is due in large part to misfits in 
the predictions associated with the Richmond-Westford 
baseline. As an example, increasing v LM from 2 x 10 2J Pa 
s to 5 x 10 23 Pas results in an increase of the reduced x 2 
by approximately 5.2 in Figure 3a. The increased misfit 
for pmlictions based on these two lower mantle viscosi- 
ties, for the Richmond-Westford baseline length rate 
alone, accounts for 50% of this increase. ( Mitrovica et 
al. [1993a] have found that far field sites, such as Rich- 
mond, have motions which are drawn more strongly to- 
ward the previously glaciated region (in this case Lau- 
rentia) as v LM is increased. This trend gives rise to 
a transition from a lengthening to a shortening of the 
Richmond-Westford baseline as v LM is increased.) 

From Figure 3 we conclude that the current VLBI 
data set of North American baseline length rates may 
be applied to assess the acceptability of specific pairings 
of ic* and Earth models. Inferences of the ice cover over 
North America, and the rheology below the region, have 
been a source^of active interest and contention in the 
geophysical literature (see below). In the following we 
consider a number of ice model/Earth mode] combina- 
tions which have been proposed for the North American 
region on the basis of a variety of independent geophys- 
ical and geological constraints. The defining character- 
istics of four particular models are given in Table 2. 

The column labelled T&P in Table 2 summarizes the 
model of Tushingkam and Peltier [1991, 1992] which 
combines the standard Earth model and the ICE-3G 


NR AO 140'- Westfofd 


Richmond - Wes I ford 



Fig. 4. The predicted baseline length rates for six base- 
lines (as labelled) computed using the ice loading history 
adapted from the ICE-3G deglaciation chronology. The 
solid, dotted and dashed lines are generated by varying, 
respectively, v LM , v UM and LT of the standard model. 
The shaded region on each frame represents the VLBI 
determined constraint (±lcr) on the baseline length rate. 


deglaciation chronology; they have argued that this com- 
bination satisfies a globally distributed data base of 
RSL curves, including those collected in North Amer- 
ica. The second model is that of Caihles [1975] (CAT) 
and includes an ice history that is similar to the ICE-1 
model of Peltier and Andrews [1976]. The CAT model 


TABLE 2. Summary of Models* and Results 
Parameter Model Values 


T&P CAT N&L DMVC Null 


VUM*! 
LT (km) 

1.0 

2.0 

120.0 

1.0 

0.9 

70.0 

0.5 

30.0 

100.0 

1.0 

1.0-10.0 

120.0 


Baseline 


Reduced ; 

K 2 


N2-W 

0.3 

1.7 

0.1 

0.1 

0.1 

R-W 

0.4 

0.0 

8.4 

0.1 

0.1 

G-R 

0.2 

0.1 

0.1 

0.1 

0.0 

G-W 

0.2 

0.5 

1.7 

0.0 

0.0 

A-G 

0.4 

3.3 

0.1 

0.2 

0.9 

C-Nl 

0.1 

0.1 

0.0 

0.0 

0.0 

All 19 

3.8 

7.9 

13.7 

2.4 

3.3 


*Ice models are discussed in text. 

** vum and i 'lm are m units of 10 21 Pa s. 
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was meant to satisfy the constraints imposed by the 
sea level record in North America, and the decay time 
inferred from that record. In contrast to these infer- 
ences, Nakada and Lambeck [1991] (N&L) have argued 
for a large (factor of 60) viscosity jump between the 
upper and lower mantle. The N&L model was based 
on fits to sea level curves in the Hudson Bay region, 
and their global ice model (ARC3 4- ANT3) adopts the 
ICE-1 chronology over North America. The final model, 
labelled DMVC (Deep Mantle Viscosity Contrast), de- 
notes a model characterized by an order of magnitude 
increase in viscosity at 1900 km depth in the mantle. 
Above 1900 km depth the viscosity is 10 21 Pa s. This 
viscosity model was proposed by Mitrovica and Peltier 
[1992] as one which might reconcile the sea level record 
over Hudson Bay, and the independent inferences of 
viscosity provided by data associated with mantle con- 
vection (see Mitrovica and Peltier , 1992, for details). 
The Mitrovica and Peltier [1992] calculations adopted 
the ICE-3G deglaciation history. 

The remaining rows of Table 2 provide the total re- 
duced x 7 f° r the VLBI determined length rates in Table 
1 relative to predictions based on the models described 
above (see row “All 19”), as well as the contribution 
to this reduced from predictions for each of the six 
baselines considered in Figure 4. The lowest (total) 
reduced x 7 values are obtained for the models DMVC 
and T&P. These reduced x 3 values are, however, in- 
distinguishable, at the 95% confidence limit, from each 
other and from the value for the null hypothesis. 

As a consequence of its low reduced x 7 the combi- 
nation of ice history and Earth rheology proposed by 
Tushingham and Peltier [1991,1992] cannot be ruled 
out by the VLBI constraints. The pairing proposed 
by Cathles [1975] yields a reduced \ 7 about double the 
T&P value; the results in Table 2 suggest that the CAT 
model performs particularly poorly in regard to predic- 
tions of length rates for the NRAO 140 ’ - Westford and 
Algonquin Park - Gilmore Creek baselines. Figures 3 
and 4 indicate that the reduced \ 2 residual might be 
lowered if the thickness of the lithosphere in the Cath- 
les [1975] model was increased. Indeed, an increase of 
LT to 120 km in the CAT model has been found to 
lower the reduced y 3 value below 4. 

The total reduced x 7 computed for the N&L ice his- 
tory/Earth rheology is relatively much higher. The pre- 
dicted length rate of the Richmond - Westford baseline 
is a primary contributor to the misfit associated with 
this model. Indeed, the Nakada and Lambeck [1991] 
model predicts a shortening of approximately 2 mm/yr 
for the Richmond - Westford baseline length rate, which 
differs significantly from the observational constraint. 
The location of Richmond, in the far field of the Lau- 
rentide ice sheet, suggests that the viscosity model pro- 
posed by these authors would almost certainly have to 
be altered in order to satisfy the VLBI determined con- 
straint on the Richmond - Westford baseline. 

The reduced \ 2 associated with the model DMVC is 
about 25% lower than the value for the null hypothe- 
sis. The result indicates that a relatively large viscosity 
contrast with depth within the lower mantle is not ruled 
out by the VLBI determined constraints on the baseline 
length rates listed in Table 1. Furthermore, it indicates 
that the sensitivity of the reduced y 3 , to variations in 
i/ LM (Figure 3), is more accurately identified as a sensi- 
tivity to variations in viscosity above 1900 km depth in 
the lower mantle (see also Mitrovica et al., 1993a). 
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